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We perform axisymmetric simulations for gravitational collapse of a massive iron core to a black 
hole in full general relativity. The iron cores are modeled by F = 4/3 equilibrium polytrope for 
simplicity. The hydrodynamic equations are solved using a high-resolution shock capturing scheme 
with a parametric equation of state. The Cartoon method is adopted for solving the Einstein 
equations. Simulations are performed for a wide variety of initial conditions changing the mass (« 
2.O-3.OM0), the angular momentum, the rotational velocity profile of the core, and the parameters 
of the equations of state which are chosen so that the maximum mass of the cold spherical polytrope 
j^-^ , is « 1.6Mq. Then, the criterion for the prompt black hole formation is clarified in terms of the 

' mass and the angular momentum for several rotational velocity profile of the core and equations of 

state. It is found that (i) with the increase of the thermal energy generated by shocks, the threshold 
mass for the prompt black hole formation is increased by 20-40%, (ii) the rotational centrifugal 
s_j ' force increases the threshold mass by < 25%, (iii) with the increase of the degree of differential 

rotation, the threshold mass is also increased, and (iv) the amplification factors shown in the results 
■ (i)-(iii) depend sensitively on the equation of state. We also find that the collapse dynamics and the 

structure of the shock formed at the bounce depend strongly on the stiffness of the adopted equation 
of state. In particular, as a new feature, a strong bipolar explosion is observed for the collapse of 
rapidly rotating iron cores with an equation of state which is stiff in subnuclear density and soft in 
supranuclear density. Gravitational waves are computed in terms of a quadrupole formula. It is also 
. found that the waveform depends sensitively on the equations of state. 
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. I. INTRODUCTION 

O ; 

The black hole is one of the most striking and intriguing objects predicted in general relativity. A wide variety of 
recent observations have shown that black holes actually exist in the universe [1]. Among several types of the black 
holes, the existence of the stellar-mass black holes has been tightly confirmed. So far, about 20 stellar-mass black 
holes for which the mass is determined within a fairly small error have been observed in binary systems of our Galaxy 
\ and the Large Magelanic Clouds [2]. Such black holes are believed to be formed from stellar core collapse of massive 
stars. This fact stimulates the theoretical study for clarifying the physics of gravitational collapse and formation 
mechanism of the stellar-mass black holes. The formation of the black hole through the gravitational collapse is a 
• >"H , highly nonlinear and dynamical phenomena. Therefore, numerical simulation in full general relativity is the unique 
y\ ' approach to this problem. 

Realistic simulations for the formation of stellar-mass black holes in massive rotating stellar core collapse is also 
increasingly important due to its possible association with gamma-ray bursts [3] . Recent observations have indicated 
that at least long-duration gamma-ray bursts are of cosmological origin [4] and associated with rotating stellar core 
collapse [5,6], probably to a black hole surrounded by a massive disk as suggested by [7-9]. Recent numerical analyses 
have also shown that if a progenitor of the collapse is massive and the angular momentum is large enough, a black 
hole surrounded by a massive disk will be indeed formed [10-12]. However detailed simulations have not yet been 
done. 

Stellar-mass black holes can be formed through the gravitational collapse of a degenerate iron core in excess of the 
Chandrasekhar mass [13,14]. It is well known that stars whose initial mass are larger than « 10M Q evolve to form a 
core mainly composed of iron group elements [14]. Since the iron is the most stable nuclei and it does not generate 
energy by nuclear burning, the iron core contracts gradually. Accordingly, the central temperature, T Cl and central 
density, p c , rise to be T c > 10 10 K and p c > 10 9 g/cm 3 , resulting in the photo-dissociation of iron to lighter elements 
and subsequent electron capture that reduce the entropy of electrons. As a result, the adiabatic index T s decreases 
below 4/3 and the iron core is destabilized to collapse. If the mass of the iron core is much larger than the maximum 
neutron star mass, a black hole will be formed soon after the collapse. However, the threshold mass of the iron core 
for prompt formation of a black hole has not been clarified yet. Note that due to the contribution of the thermal 
pressure and rotation, the iron core can be much larger than the maximum allowed neutron star mass of ~ 2M® for 
a sufficiently massive and rotating star. 



1 



Since massive stars in nature are rapidly rotating in general [15], it is necessary to explore the gravitational collapse 
of a rotating star in full general relativity for clarifying the black hole formation. Since Nakamura [16] and his 
collaborators [17] first presented a series of preliminary numerical simulations, a number of simulations for rotating 
stellar collapse to a black hole have been performed in full general relativity [18-23]. However, the initial conditions 
and equations of state in the previous works are not very realistic for modeling the stellar core collapse. Thus, the 
simulations with a realistic setting remain to be an unsolved issue in general relativity. 

Realistic simulations of rotating stellar core collapse have been performed intensively in the framework of Newtonian 
gravity [24-33] . Most of these studies mainly aim at clarifying the effect of the rotation on the dynamics of neutron 
star formation and gravitational waveforms from it. In particular, a comprehensive parameter study sweeping through 
various values of the stiffness of a parametric equation of state as well as rotational parameters was performed in [29] . 
It was shown that the dynamics of the collapse and resulting gravitational waveforms depend strongly not only on 
rotation but also on the stiffness of equations of state. Dimmelmeier et al. [34] extended the aforementioned study 
to general relativistic case using a conformal flatness formalism [35]. Fully general rclativistic numerical studies of 
neutron star formation have been recently performed [36,37]. As shown in [34,36], the general relativistic effects 
significantly modify the dynamics of the collapse even in the formation of neutron stars. 

Taking into account the present status described above, in this paper, we study a criterion for prompt black hole 
formation in the iron core collapse performing fully general relativistic simulations. The iron cores are modeled by 
T = 4/3 polytropes in equilibrium for simplicity. The major purpose of this paper is to clarify the threshold mass 
of the iron core for the prompt black hole formation and its dependence on the angular momentum, the rotational 
velocity profiles of the iron core, and the equations of state. To clarify the dependence on the equations of state in a 
clear manner, we adopt a parametric equation of state following previous papers [36,38]. 

The simulations are performed assuming that the collapse proceeds in an axisymmetric manner. This assumption 
is reasonable as far as the progenitor of the collapse is not very rapidly and highly differentially rotating (e.g., [38]). 
In this paper, we do not adopt such progenitor that are likely to be dynamical unstable against nonaxisymmetric 
deformation (cf. Sec. IV D 2 for discussion). Although for several models, the rotational kinetic energy is so large 
that the outcome formed in the collapse may be secularly unstable against nonaxisymmetric deformation, the secular 
time scale is much longer than the dynamical time scale of the core collapse. Hence, the collapse will proceed in 
an approximately axisymmetric manner in the time scale of interest. On the other hand, rapidly and differentially 
rotating stellar core collapse has to be studied in the three-dimensional simulation. Such simulation was recently 
performed and the detailed results are shown in a companion paper [38]. 

This paper is organized as follows. In Sec. II, we briefly review our formulation for general relativistic hydrodynamic 
simulations, equations of state, and a quadrupole formula adopted in the present paper. In Sec. Ill, we describe the 
initial conditions. A detail of computational setting is also described. Sec. IV presents the numerical results, 
emphasizing the threshold mass for the prompt black hole formation and its dependence on the angular momentum 
and the adopted equations of state. Gravitational waveforms emitted in the neutron star formation are also shown. 
Sec. V is devoted to a summary. Throughout the paper, we adopt the geometrical units G = c = 1 where G and c are 
the gravitational constant and speed of light, respectively. The Latin indices i, j, k, ■ ■ ■ denote the spatial components 
of x, y, and z, and the Greek indices \i ■ ■ ■ denote the spacetime components. 

II. NUMERICAL IMPLEMENTATION 
A. Brief summary of formulation and numerical method 

We perform fully general relativistic simulations for rotating stellar core collapse in axial symmetry using the same 
formulation and numerical method as those presented in [39], to which the reader may refer for details of basic 
equations and successful test simulations. 

In the 3+1 formulation, the metric can be written in the form 

ds 2 = {-a 2 + f3 k f3 k )dt 2 + 2/3 k dtdx k + -y ij dx i dx j , (1) 

where a, f3 k , and 7^ are the lapse function, the shift vector and metric in 3D spatial hypersurface, respectively. The 
extrinsic curvature is defined by 

(d t - fy^ij = -2aK ijt (2) 
where £p is the Lie derivative with respect to (3 k . 
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As in the series of our papers, we evolve <j> = log(det7y)/12, 7^- = e^^jij, Aij = e~ A< ^(Kij — "fijK k k ), and 
trace of the extrinsic curvature K k together with three auxiliary functions Fi = S^ k dj^fik with an unconstrained 
free evolution code as in [40,20,41-44,39]. The Einstein equations are solved in Cartesian coordinates. To impose 
axisymmetric boundary conditions, the Cartoon method [45] is used with the grid size N x 3 x N in (x, y, z) which 
covers a computational domain as < x < L, < z < L, and —A < y < A. Here, TV and L are constants and 
A = L/N. 

The fundamental variables for the hydrodynamics are p : rest mass density, e : specific internal energy, P : pressure, 
: four velocity, and 

v' = ¥- = £ t . (3) 

As the variables to be evolved in the numerical simulations, we define a weighted density, /?*, a weighted four- velocity 
u^, and a specific energy density as 

= puie 6 *, 
Ui = hm, 

P 

e = hw , (4) 

pw 

where w = cm*, and h = 1 + e + P/ p. From these variables, the total baryon rest mass and angular momentum of 
the system, which are conserved quantities in an axisymmetric spacetime, can be defined as 



f* = J (fxp^, (5) 
J = J d z xp*u v . (6) 



The general relativistic hydrodynamic equations are solved using a so-called high-resolution shock-capturing scheme 
[46] on the y = plane with the cylindrical coordinates (x, z) (in Cartesian coordinates with y = 0). Details about 
our numerical scheme are described in [39] 

We neglect effects of viscosity and magnetic fields. The time scale of dissipation and angular momentum transport 
due to these effects are much longer than the time scale of collapse ~ 100 ms, unless the magnitude of viscosity or 
magnetic fields is extremely large [47]. Thus neglecting them is an appropriate assumption. 

As the slicing condition we impose an "approximate" maximal slicing condition in which K k k w is required [40] . 
As the spatial gauge, we adopt a dynamical gauge condition [48] in which the equation for the shift vector is written 
as 

d t (3 k =^ kl (F l + Atd t F l ). (7) 

Here, At denotes the time step in numerical computation [21]. Note that in this gauge condition, (3 l obeys a hyperbolic- 
type equation for a sufficiently small value of At because the right-hand side of the evolution equation for Fi contains 
a vector Laplacian term [43] . It has already been illustrated that stable simulations for rotating stellar collapse and 
merger of binary neutron stars are feasible in this gauge [21,42]. 

An outgoing- wave boundary condition is imposed for hij(= jij — Sij), Aij, and Fi at the outer boundaries of the 
computational domain. The condition adopted is the same as that described in [43]. On the other hand, for <j> and 
K k k , other types of outer boundary conditions are imposed as r (j> — const and K k = 0, respectively. 

Existence of a black hole is determined using an apparent horizon finder developed in [49] . We compute the apparent 
horizon mass Mah which is defined as [50] 



Mah = V^ (8) 

where A denotes area of an apparent horizon. 

During the numerical simulations, conservation of the Arnowitt-Deser-Misner (ADM) mass, Madm, and the angular 
momentum are monitored as code checks. The ADM mass is defined by 



M ADM = J d 3 ^{phw 2 - P) + ^- (l^ - 2 -{K k )^j - 



**** 



(9) 
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where ip = . 

For the analysis of numerical results, we define a rest-mass distribution m*(j) [19], which is the integrated baryon 
rest mass of fluid elements with the specific angular momentum less than a given value of j = u v : 



m*(j) = 27r/ p^r 2 drd{cos9). (10) 
Jj'<j 

Similarly, a specific angular momentum distribution is defined according to 

J(j)=2n [ p*j'r 2 drd(cos6). (11) 
Jj'<i 

These distribution functions are preserved in axisymmetric spacetimes of ideal fluid. Gauge independence and preser- 
vation of these distribution functions in axial symmetry can be proven by the hydrodynamical equations 

«M + ^W)=, (13) 

where the index I denotes the component of w and z. 

From these distribution functions, we define a spin parameter distribution as 

«•«> s <»» 

This may be approximately regarded as a spin parameter of the inner region of the core composed of fluid elements 
with the specific angular momentum less than j. 



B. Equations of state 



Model £i £2 pnuc(g/cm 3 ) Tth M max 

a L32 2^25 2.0 x 10 14 L32 1.623 

b 1.30 2.5 2.0 x 10 14 1.30 1.600 

c 1.30 2.22 1.0 x 10 14 1.30 1.599 

d 1.28 2.75 2.0 x 10 14 1.28 1.597 



TABLE I. Adopted sets of (ri, T2, p n uc, Tth). The values of Ti, T2, and p nuc are chosen so that the maximum ADM mass of 
a cold spherical polytrope for each set becomes w 1.6M . 

During dynamical evolution, a parametric equation of state is adopted following Miiller and his collaborators 
[51,29,34]. In this equation of state, one assumes that the pressure consists of the sum of polytropic and thermal parts 
as 

P = P P +P th . (15) 

Pp denotes the cold (zero temperature) nuclear equation of state and is given by P P = K P (p)p r ^ where K P and T 
are functions of p. In this paper, we follow [34] for the choice of Kp(p) and T(p): For the density smaller than the 
nuclear density, p n uc, we set T = Ti < 4/3, and for p > p nuc , T — T 2 > 2. Namely, 

J Kip T \ p < p nuc , . . 

Pp -\K 2 p^, p> Pnuc , ( 16 ) 

where K\ and K2 are constants. Since Pp should be continuous at p = p nuc , the relation, K2 — Kip^~ T2 , is required. 
Following [29,34], the value of K\ is fixed to be 5 x 10 14 cgs. With this choice, the polytropic part of the equation of 
state for p < p nuc , in which the degenerate pressure of electrons is dominant, is approximated well. Since the specific 
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FIG. 1. The pressure, P, as a function of the density, p, for cold equations of state with the parameters listed in Table I. 



internal energy should be also continuous at p = p nuc , the polytropic part of the specific internal energy, e P , is denned 

as 



£p 



Ki r i 

r<° 1 , P < P» 



r — i 

l r 2 - i + (r!-i)(r 2 -i) ' p - p ™- 



With this setting, a realistic equation of state for high-density, cold nuclear matter is mimicked. 

For realistic simulations of stellar core collapse, it would be better to adopt realistic equations of state [52,53]. 
However, in the realistic equations of state, many micro-physical processes are simultaneously taken into account 
together. In such cases, it is not easy to extract an important element responsible for an output of numerical 
simulations. With the parametric equations of state, on the other hand, one can systematically investigate dependence 
of the dynamics of stellar core collapse on the equations of state by changing their own parameters. Therefore, as a 
first step to the more realistic simulations of rotating stellar core collapse, we adopt the parametric equations of state. 
We now plan to perform simulations with realistic equations of state. Some discussion comparing the parametric 
equations of state and a realistic equation of state is presented in appendix A. 

There is four parameters in the parametric equations of state, namely (T\, r 2 , Pnuc, Tth)- In this paper, we choose sets 
of (Ti, r 2 , Pnuc) so that the maximum allowed ADM mass of the cold spherical polytrope becomes an approximately 
identical value as MADM,max ~ 1.6M©. Note that this value is larger than the mass of neutron stars in binary 
neutron stars accurately determined, and thus, a reasonable choice [54]. Following [29,34], we typically set p nuc = 
2.0 x 10 14 g/cm 3 (pi4 = pnuc/(10 14 g/cm 3 ) = 2.0). For Ti, we choose the three values of 1.32, 1.30, and 1.28. Requiring 
that the maximum allowed ADM mass should be w 1.6M Q , the values of T 2 are determined to be 2.25, 2.50, and 
2.75, respectively. To investigate the dependence of the output physics on the value of p nuc , we also pick up a case 
with p nuc = 1.0 x 10 14 g/cm 3 (p 14 = 1.0) and set (Fi,r 2 ) = (1.30,2.22) for comparison. 

To summarize, the parameter sets adopted in this paper are (ri,r 2 ,pi 4 ) = (1.32,2.25,2.0), (1.30,2.50,2.0), 
(1.30,2.22,1.0), and (1.28,2.75,2.0) which are referred to as equations of state 'a', 'b', V, and 'd', respectively 
(cf. Table I) . Although the maximum allowed mass of the cold spherical neutron star in equilibrium is approximately 
identical for all of the equations of state, the difference in the values of (ri, T 2 , p nuc ) yields a significant variation in 
the collapse dynamics and in the criterion for prompt black hole formation. 

Figure 1 shows relations between the pressure and the density for each set of (ri, T 2 , p nuc ). It is found that for the 
smaller value of Ti, the depletion fraction of the pressure for p < p nuc is increased by a large factor. It is also worthy 
to note that the equation of state 'c' is stiffcr than 'b' in the density range between 10 and «2x 10 15 g/cm 3 . These 
result in significant difference in the collapse dynamics and the threshold mass for the prompt black hole formation. 

Pth is related to the thermal energy density, e t h = £ — £p , as 

Pth = (r th - l)pe th . (18) 
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In this paper we set T t h = Ti for simplicity. As shown in [36], the collapse dynamics depends rather weakly on I\h 
as far as it is in the range 4/3 < T t h ^ 5/3, besides the fact that the oscillation amplitude of a formed protoneutron 
star depends slightly on it. 

To prepare rotating stars in equilibrium as initial conditions, we use the polytropic equations of state with T = 4/3 

P = K p 4 / 3 = (K dcg + K T )p 4 '\ (19) 

where following [34], K<± e g is set to be 5 x 10 14 cm 3 /s 2 /gr 1 / 3 , with which a soft equation of state governed by the 
electron degenerate pressure is approximated well [13]. Here, K& eg and K\ are related by K\ ~ K^p^ 3 ~ Tl where 
we set po = 1 g/cm 3 . The extra pressure of Pt = K^p 4 / 3 denotes the pressure generated by the ideal gas pressure 
and the radiation pressure that are nonzero only for the finite temperature. The values of K and Kt adopted in the 
present paper are described in Sec. III. To induce the collapse, we slightly decrease the value of the adiabatic index 
from T = 4/3 to T\ < 4/3 at t = 0. This implies that at t = 0, Pp = Kd CE p Tl and P t h = KtP Ti , respectively. 



C. Quadrupole formula 

In the present work, gravitational waveforms are computed using a quadrupole formula described in [55,36]. In 
quadrupole formulas, only the +-modc of gravitational waves with I = 2 is nonzero in axisymmetric spacetime and it 
is written as 

^quad = J«(trot)-4s(frot) ^2 ^ (2Q) 

where Iij denotes a quadrupole moment, 7y its second time derivative, and t re t a retarded time. 

In fully general relativistic and dynamical spacetime, there is no unique definition for the quadrupole moment and 
nor is for 1^ . Following a previous paper [55] , we choose the simplest definition as 

hj = j p*x l x 3 d 3 x. (21) 

Then, using the continuity equation of the form 

dtp* + d t (p*v l ) =0, (22) 

the first time derivative can be written as 



Iij 



= J p^x 1 +x l v 3 )d 3 x. (23) 



To compute 1^, the finite differencing of the numerical result for 1^ is carried out. 
In the following, we present 

^(irot) = Izz(trct) - Ixx(trct), (24) 

in the quadrupole formula. This provides the amplitude of I = 2 mode measured by an observer located in the most 
optimistic direction (i.e., in the equatorial plane). 

The energy power spectrum of gravitational waves is given in [56] as 

d -^ = \i^)f(\h + (m (/>0), (25) 

where h + {f) denotes the Fourier transform 

h+(f) = J ^ uad (t)e 2m/t ^, (26) 

and the bracket denotes the angle-averaged value. This can be expressed in the present case as (e.g. [57]) 

<IM/)| 2 } = ~|^(/)f, (27) 
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where 



Mf) = J A 2 {t)e 2 ^ l dt. (28) 
Thus, the effective amplitude of gravitational waves observed in the most optimistic direction is denoted by 

As indicated in [55] , it is possible to compute gravitational waves from oscillating and rapidly rotating neutron stars 
of high values of compactness fairly accurately with the present choice of I{j , besides possible systematic errors for the 
amplitude of order M/R. For the stellar core collapse in which the outcomes are protoneutron stars of M/R ~ 0.1-0.2, 
it is likely that the wave amplitude is computed within an error of ~ 10-20%. The phase of gravitational waves will 
be computed very accurately as indicated in [55]. For the stellar core collapse to a black hole, on the other hand, 
the quadrupolc formula will be no longer valid because the value of M/R is high and, moreover, quasi-normal mode 
ringing of the formed black hole is not taken into account. 

III. INITIAL CONDITIONS AND COMPUTATIONAL SETTING 

A. Initial conditions 

Recent numerical study [58,59] for stellar evolution of very massive and low metallicity stars from main-sequence to 
the pre- iron-core-collapse suggests that initially massive stars evolve to form an iron core of mass > 2Mq-3M q . Taking 
these fact into account, we consider a wide range of the mass for the progenitor of the collapse as 2M < M < 3M Q . 

It should be addressed that such large mass of the iron core may not be a special product of low metallicity or large 
progenitor mass. Recently, Hirschi ct al. [60] study presupernova evolution of rotating massive star of solar metallicity 
and show that mass of the iron core at the onset of the collapse depends strongly on the treatment of both convection 
and rotation. They find that the rotation significantly increases (by a factor of ~ 1.5) the core mass: An iron core of 
mass of w 2.0Mq is formed from the ZAMS (zero-age main sequence ) of 15M by including the rotation while its 
non-rotating counterpart yields an iron core of mass w 1.5M©. This result is quite different from that in [61]. This 
suggests that mass of the iron core at the onset of collapse may depend sensitively on the chemical abundance, and 
the stellar rotation as well as the magnitude of the viscosity and the mixing length which are not well understood. 

The central density and the central temperature of such very massive iron cores are p c > 5 x 10 9 g/cm 3 and 
T c ~ 10 10 K, respectively [59]. Thus, we set the central density of the initial conditions to be 10 10 g/cm 3 . The 
electrons under such a high density are extremely degenerate even with T w 10 10 K, since ratio of the Fermi energy 
of free electrons, ep, to thermal energy, fc^T, is much larger than unity [13] according to 

Sf" 10 ^) (lOio g/cm 3 ) (iO^k) ' (30) 

Therefore, the main contribution of the pressure comes from the electron degenerate pressure that is denoted by the 
polytropic form as Pdcg = -^dcgP 4 ^ 3 - 

In addition to the electron degenerate pressure, the thermal and radiation pressure should be taken into account. 
The adiabatic index relevant for them in the iron core may be close to 4/3 at the onset of the collapse since the 
photo-dissociation of the irons are accelerated with increasing the density to reduce the pressure. Hence, we put the 
gas and radiation pressure together into the polytropic form as 

P T = -Pgas + -Prad = K T p$ . (31) 

With these assumptions, rotating polytropes in equilibrium with the polytropic index V = 4/3 and the polytropic 
constant K = i^dog + Kt [see Eq. (19)] are given as the initial models of rotating iron cores. 

The iron core model adopted in this paper is nothing more than a simplified one. The nature is more complicated 
and the iron core of the evolved massive star are not simple polytropes [33]. We will nevertheless use the term 'iron 
core' in the following, to emphasizing that we aim at clarifying the threshold mass of the iron core for the prompt 
black hole formation. 
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105/120 


2.449 


2065 


f on , , 1 n — 3 

5.20 x 10 


0.792 


3.20 


D3 


^7 i~\ * , 1 14 

7.0 x 10 


100/120 


2.459 


2162 


6.49 x 10 


0.888 


3.55 


Pi A 


/ .U X 1U 


nn /i on 
yu/ izu 


O A 70 
Z.4 / z 


Zoo / 


o.ZO X 1U 


l.UUO 


Q Q7 

o.y / 


Do 


/.0 X 10 


on / 1 on 
8U/ IzU 


2.4/0 


20/ 9 


8.80 X 10 


1 n a k 
1.040 


a 1 n 
4.10 




/.ZO X 1U 


1 1 n /1 on 
1 iu / izu 


z.OOo 


zuiy 


O.00 X 1U 


U.OOz 


7n 
z. /U 


E25 


7.25 X 10 


1 n cr / 1 on 

105/120 


2.580 


2103 


5.20 x 10 


0.780 


3.19 


E3 


7.25 x 10 14 


100/120 


2.590 


2198 


6.49 x 10~ 3 


0.873 


3.55 


E4 


7.25 x 10 14 


90/120 


2.604 


2428 


8.25 x 10~ 3 


0.989 


3.97 


E5 


7.25 x 10 14 


80/120 


2.608 


2726 


8.86 x 10 -3 


1.03 


4.09 


FO 


7.5 x 10 14 


1 


2.672 


1909 


0.0 


0.0 


0.0 


Fl 


7.5 x 10 14 


117/120 


2.680 


1951 


1.17 x 10~ 3 


0.361 


1.55 


F2 


7.5 x 10 14 


110/120 


2.700 


2056 


3.66 x 10 -3 


0.641 


2.71 


F3 


7.5 x 10 14 


100/120 


2.723 


2238 


6.49 x 10~ 3 


0.858 


3.55 


F4 


7.5 x 10 14 


90/120 


2.738 


2471 


8.26 x 10 -3 


0.972 


3.97 


F5 


7.5 x 10 14 


80/120 


2.742 


2773 


8.86 x 10~ 3 


1.010 


4.10 


GO 


7.75 x 10 14 


90/120 


2.874 


2510 


8.24 x 10~ 3 


0.957 


3.96 


G5 


7.75 x 10 14 


80/120 


2.879 


2818 


8.85 x 10 -3 


0.993 


4.09 


HO 


8.0 x 10 i4 


1 


2.940 


1972 


0.0 


0.0 


0.0 


H5 


8.0 x 10 14 


80/120 


3.016 


2864 


8.84 x 10~ 3 


0.978 


4.10 



TABLE II. The values of Ko, axial ratio R p /R e , ADM mass Madm, equatorial circumferential radius _R CC , ratio of the 
rotational kinetic energy, T rot , to the gravitational potential energy, W, nondimensional spin parameter q = .//Af| DM , and 
angular velocity fi a for rigidly rotating iron cores in equilibrium used as the initial condition for numerical simulation. 
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Model 


K (cgs) 


A 


Madm(M q ) 


-Rc(km) 


T rot /W 


q 


n„(i/s) 


D20dl 


7.0 x 10 14 


1.0 


2.437 


1912 


3.54 x 10 -3 


0.647 


3.06 


D22dl 


7.0 x 10 14 


1.0 


2.441 


1924 


4.12 x 10" 3 


0.698 


3.29 


D23dl 


7.0 x 10 14 


1.0 


2.445 


1936 


4.69 x 10~ 3 


0.746 


3.50 


D25dl 


7.0 x 10 14 


1.0 


2.450 


1950 


5.26 x 10~ 3 


0.791 


3.69 


D15d05 


7.0 x 10 14 


0.5 


2.424 


1858 


1.75 x 10~ 3 


0.440 


2.84 


D17d05 


7.0 x 10 14 


0.5 


2.430 


1866 


2.63 x 10~ 3 


0.541 


2.84 


D20d05 


7.0 x 10 14 


0.5 


2.437 


1874 


3.52 x 10" 3 


0.626 


3.47 


D25d05 


7.0 x 10 14 


0.5 


2.444 


1881 


4.42 x 10 -3 


0.702 


4.21 


H5dl 


8.0 x 10 14 


1.0 


3.019 


2190 


9.03 x 10~ 3 


0.978 


4.75 


H5d05 


8.0 x 10 14 


0.5 


3.020 


2054 


8.98 x 10~ 3 


0.941 


6.31 



TABLE III. The values of K , A = zu d /R e , Madm, Re, T rot /W, q = J/M^ DM , an d angular velocity at the rotational axis 
fi a for differentially rotating initial models. 



We adopt the values of K u = K Q /10 1A (cgs) to be 6.0, 6.5, 6.75, 7.0, 7.25, 7.5, 7.75, and 8.0, which we refer to 
models A, B, C, D, E, F, G, and H, respectively (cf. Table II). In the Newtonian case, mass of the T = 4/3 spherical 
polytrope is given by M w 4.555(X /G) 3/2 . Thus, for the adopted values of K , the ADM mass of the initial spherical 
iron cores are Mabm/Mq w 2.0, 2.2, 2.3 2.4, 2.6, 2.7, 2.9, and 3.0, respectively (see Table II). These values are larger 
than the maximum allowed ADM mass, w 1.6M©, of the spherical cold polytrope adopted in this paper. Namely, in 
the absence of shock heating and rotation, the core collapses to a black hole. 

For our choice of the iron core mass, the corresponding helium core mass is M^/Mq « 17, 22, 25, 30, 32, 35, 40, 
and 42 for models A, B, C, D, E, F, G, and H, respectively, according to the calculation by Umeda and Nomoto [59]. 
Note that these values are larger than a critical value of helium core mass for direct black hole formation w 15M 
estimated in [62,63]. 

The velocity profiles of equilibrium rotating cores are given according to a popular relation [64,65] 

u*u v = ru 2 d (n a - ft), (32) 

where Q a denotes the angular velocity along the rotational axis, and m d is a constant. In the Newtonian limit, the 
rotational profile is written as 

2 

n = n a 2 Wd 2 . (33) 

Thus, vud controls the steepness of differential rotation. In this paper, we pick up mainly the rigidly rotating models 
in which Wd — > oo. For illustration of the effect of differential rotation, we select differentially rotating models with 
A = Wd/ R e = 1 and 0.5 (see table III). 

For rigidly rotating initial models, we choose the axial ratio R p /R e of polar radius, R p , to equatorial radius, R e , as 1 
(spherical configuration), 119/120, 115/120, 110/120, 105/120, 100/120, 90/120, and 80/120. The models with these 
axial ratios are referred to as models X0, XI, X15, X2, X25, X3, X4, and X5, where X denotes A-H. For example, a 
model with K 1A = 7.0 and R p /R e = 80/120 is abbreviated as D5. Note that for models with R p /R e = 80/120, the 
angular velocity at the equatorial stellar surface is nearly equal to the Keplerian velocity; namely, a rapidly rotating 
initial condition near the mass shedding limit is chosen for this case. On the other hand, differentially rotating initial 
models are chosen only for a selected set of parameters with K u = 7.0 and 8.0, since the purpose in this paper is 
to clarify the effect of differential rotation, comparing the results with those for the rigidly rotating initial models of 
nearly identical values of mass and angular momentum. 

In Tables II and III, several fundamental quantities for the models adopted in the present numerical computation 
are listed. Here, q is a nondimensional spin parameter defined by J/M^ DM . T Iot and W are the rotational kinetic 
energy and the gravitational potential energy, and defined according to [66,67] 



T rot = t I cPxp.u v n, (34) 



j3. 



W= M* + J p^ed'x - Madm + T rot , (35) 
where W is defined to be positive. Note that for the rigidly rotating case, the maximum value of T TOt /W is « 0.009. 
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For the rigidly or weakly differentially rotating stars, the stability against nonaxisymmetric perturbation is deter- 
mined by T TOt /W [68]. A star with T rot /W > 0.27 and T Iot /W > 0.14 will become unstable against the nonaxisym- 
metric dynamical and secular instabilities, respectively. Here, the secular timescale is much longer than the collapse 
timescale, so that the dynamical instabilities are only relevant during collapse and bounces. Since we assume that the 
rotating iron core collapse proceeds in an axisymmetric manner, it is important to ensure that the iron core should 
not spin up to be dynamically unstable against nonaxisymmetric deformation. We will discuss the spin-up and the 
stabilities of the inner cores formed after the bounce in Sec. IV D 2. 

For the differentially rotating case with a small value of A, it is possible to make equilibrium states with T rot /W 
0.009. With such an initial condition, the collapsing core may form a differentially rotating object of a highly 
nonspherical shape and of a high value of T rot /W [38,30]. It is also known that rapidly rotating neutron stars of a high 
degree of differential rotation is dynamically unstable against nonaxisymmetric deformation even for T rot /W = 0(0.01) 
(e.g., [69] and references therein). To follow such collapse, a nonaxisymmetric simulation will be necessary. In this 
paper, we do not choose such initial conditions and focus only on dynamically stable cases against nonaxisymmetric 
deformation. The collapse of highly differentially rotating initial conditions is studied in three-dimensional simulations 
in [38]. 

B. Computational Settings 

The central density increases from 10 10 g/cm 3 to > 10 15 g/cm 3 during the collapse. This implies that the char- 
acteristic length scale of the system varies by a factor of ~ 100. To compute such a collapse accurately saving the 
CPU time efficiently, a regridding technique as described in [19,36] is helpful. The regridding is carried out whenever 
the characteristic radius of the collapsing star decreases by a factor of a few. At each regridding, the grid spacing is 
decreased by a factor of 2. All the quantities in the new grid are calculated using the cubic interpolation. To avoid 
discarding the matter in the outer region, we also increase the grid number at the regridding. 

For the regiridding, we define a relativistic gravitational potential $ c = 1 — a c ($ c > 0). Since <fr c is approximately 
proportional to M/R, $~ 1 can be used as a measure of the characteristic length scale of the core for the regridding. 
We typically choose N at each regridding in the following manner. From t = to the time at which $ c = 0.025, we 
set N = 620 with the grid spacing Ai? e /600. At $ c = 0.025, the characteristic stellar radius becomes approximately 
one fourth of the initial value. Then, the first regridding is performed; the grid spacing is changed to the half of the 
previous one and the grid number is increased to N = 1020. Subsequently, the value of N is chosen in the following 
manner; for 0.025 < $ c < 0.05, we set N = 1020; for 0.05 < $ c < 0.1, we set N = 1700; and for 0.1 < $ c < 0.25, we 
set N — 2500, and keep this number until the termination of the simulations. For the typical cases, the physical size 
of the grid spacing is A ~ 4 km at the beginning of simulations and ~ 0.5 km at the end. 

In the case of black hole formation, <!> c approaches to 1. In this case, we carry out one more regridding at the time 
of $ c = 0.25. In this final regridding, the grid spacing is made half while keeping N = 2500. For the typical cases in 
which a black hole is formed, the physical size of the grid spacing is A ~ 4 km at the beginning of simulations and 
~ 0.25 km at the end. In this treatment, the total discarded fraction of the baryon rest-mass which is located outside 
the new regridded domains is < 4%. 

Simulations for each model with the higher grid resolution are performed for 40,000-100,000 time steps. The required 
CPU time for one model is about 30-90 hours using 8 processors of FACOM VPP 5000 at the data processing center 
of National Astronomical Observatory of Japan. 

IV. NUMERICAL RESULTS 

A. General feature of the collapse 

We first summarize the outline of the dynamics of stellar core collapse. Detailed features of the collapse will be 
presented in the subsequent subsections. In Figs. 2-5, we display evolution of the central density p c and the central 
value of the lapse function a c for selected models. As indicated in these figures, rotating stellar core collapse to a 
neutron star can be divided into three phases; the infall phase, the bounce phase, and the ringdown (or post-bounce 
oscillation) phase [26,29]. 

The infall phase sets in due to the onset of the gravitational instability of the progenitor triggered by the sudden 
softening of the equation of state which is associated with the reduction of the adiabatic index. During this phase, the 
central density, p Cl (the central value of the lapse function, a c ), monotonically increases (decreases) until it reaches 
the nuclear density, provided that the core is not very rapidly rotating initially. The inner part of the core, which 
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FIG. 2. Evolution of the central density p c (upper panel) and the central value of the lapse function a c (lower panel) for 
models D2a (long dashed curve), D3a (dashed curve), D5a (solid curve), and F5a (dotted curve). 
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FIG. 3. The same as Fig. 2 but for models Dlb (long dashed curve), D2b (dashed curve), D5b (solid curve), and F5b (dotted 
curve) . 



collapses nearly homologously with a subsonic infall velocity, constitutes the inner core. On the other hand, the outer 
region in which the infall velocity is supersonic constitutes the outer core [26,29]. 

The bounce phase sets in when the density around the central region exceeds the nuclear density. At this phase, 
the inner core decelerates rapidly due to (a) the sudden stiffening of the equation of state or (b) the strong centrifugal 
force. Hereafter, we pay attention to the case (a) since the collapse is halted due to the sudden stiffening of the 
equation of state for all the models (cf. Figs. 2-5). Because of its large inertia and large kinetic energy induced by 
the infall, the inner core overshoots its hypothetical equilibrium state. The degree of the overshooting depends on the 
mass, the amount of rotational kinetic energy and the stiffness of the equation of state. 

If its mass is not too large, the inner core experiences a bounce. The stored internal energy of the inner core at 
maximum compression is released through a strong pressure wave generated inside the inner core [26]. The pressure 
wave travels from the center to the outer region until it reaches the sonic point located at the edge of the inner core. 
Since the sound cones tilt inward beyond the sonic point, the pressure disturbance cannot travel further and forms 
a shock just inside the sonic point. During the formation of the shock, the inner core transfers its kinetic energy to 
the shock through compressional work that powers the shock [70,26]. It is important to note that the shock is formed 
at the outer edge of the inner core, and hence, the bulk of the inner core matter never undergoes shock heating and 
acceleration. Therefore, the shock heating is relatively less efficient for models with larger mass of the inner core. 

On the other hand, if the mass and inertia of the inner core at the bounce are sufficiently large, the pressure supplied 
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FIG. 4. The same as Fig. 2 but for models Die (long dashed curve), D15c (dotted dashed curve), D5c (solid curve). Results 
for models Dlb (dotted curve) and D5b (dashed curve) are also displayed for comparison. 




FIG. 5. The same as Fig. 2 but for models Did (dotted curve), D2d (long dashed curve), D3d (dotted dashed curve), and 
D5d (solid curve). 



by the sudden stiffening of the equation of state and the centrifugal force cannot halt the collapse. Then, the inner 
core will promptly collapse to a black hole without any distinct bounce (cf. models D2a, D3a, and D5a in Fig. 2 in 
which a black hole is formed). As illustrated in Sec. IV B 2, shocks do not propagate outward in such cases. 

In the ring-down phase, the inner core oscillates quasi-radially and then settles down to a quasi-stationary state, 
since the compressional work done by the inner core on the matter of the outer region leads to damping of the 
oscillation. In this phase, the amplitude of the oscillation and the strength of shocks generated by the outward 
oscillation depend on the stiffness of the equation of state for the formed protoneutron star. In the outer region, on 
the other hand, shock waves propagate and are accelerated due to the density gradient in the outer core. They sweep 
materials of the outer envelopes, and convert the infall kinetic energy into the thermal energy, which helps further 
driving the shock outward. However, if the explosion is too weak to eject sufficient matter of the progenitor star, the 
subsequent fallback of the matter into the formed protoneutron star will trigger formation of black hole. 



B. Criterion for prompt black hole formation 



In Table IV, we summarize the outcomes in the iron core collapse for all of rigidly rotating initial models listed in 
Table II and with equations of state listed in Table I. It is important to note that even for models in which a black 
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K (cgs) 



Rp/ Re 



Maum(Mq) 



BO 


6.5 x 10 


1 


2.163 


0.0 


BH 


BH 


no BH 


BH 


Bl 


6.5 x 10 14 


119/120 


2.163 


0.217 


no BH 


no BH 


no BH 


no BH 


B2 


6.5 x 10 14 


110/120 


2.185 


0.688 


no BH 


no BH 


no BH 


no BH 


"DO 

B3 


6.0 X 10 


i nn / 1 on 
1UU/ IzU 


2.203 


0.921 


no tin 


no titi 


no on 


~ ~ 1 1 1 1 
no BH 


B4 


6.5 x 10 14 


90/120 


2.215 


1.044 


no BH 


no BH 


no BH 


no BH 


B5 


6.5 x 10 14 


80/120 


2.219 


1.084 


no BH 


no BH 


no BH 


no BH 


CO 


6.75 x 10 


1 


2.286 


0.0 


BH 


BH 


BH 




CI 


6.75 x 10 14 


119/120 


2.289 


0.221 


BH 


BH 


no BH 


— 


C15 


6.75 x 10 14 


115/120 


2.298 


0.487 


BH 


BH 


no BH 




C2 


6.75 x 10 14 


110/120 


2.310 


0.675 


BH 


no BH 


no BH 




C25 


6.75 x 10 14 


105/120 


2.321 


0.807 


no BH 


no BH 


no BH 




DO 


7.U X 10 


1 


2.412 


0.0 


BH 


BH 


BH 


BH 


Dl 


7.0 x 10 14 


119/120 


2.414 


0.217 


BH 


BH 


BH 


BH 


D15 


7.0 x 10 14 


115/120 


2.425 


0.478 


BH 


BH 


BH 


BH 


D2 


7.0 x 10 14 


110/120 


2.438 


0.663 


BH 


BH 


no BH 


BH 


D25 


7.0 x 10 14 


105/120 


2.449 


0.792 


BH 


BH 


no BH 


BH 


Uo 


(.[) X 1U 


i nn 1 1 on 
1UU/ IzU 


2.4o9 


0.000 


DU 

ati 


no BH 


no BH 


no BH 


D4 


7.0 x 10 


90/120 


2.472 


1.006 


no BH 


no BH 


no BH 


no BH 


D5 


7.0 x 10 


80/120 


2.476 


1.045 


no BH 


no BH 


no BH 


no BH 


E2 


7.25 x 10 14 


110/120 


2.568 


0.652 


BH 


BH 


BH 


— 


E25 


7.25 x 10 14 


105/120 


2.580 


0.780 


BH 


BH 


no BH 


— 


E3 


7.25 x 10 


100/120 


2.590 


0.873 


BH 


BH 


no BH 






*7 oc w i nl4 
I.ZO X 10 


nn 1 1 on 

yu/izu 


o tin A 
2.D04 


a non 
0.989 


ati 


no BH 


no BH 




E5 


7.25 x 10 


80/120 


2.608 


1.03 


BH 


no BH 


no BH 




r U 


/.0 X 1U 


i 
l 


2.0 1 1 


U.U 




JdH 


BH 


BH 


Fl 


7.5 x 10 


119/120 


2.680 


0.361 


BH 


BH 


BH 


BH 


F2 


7.5 x 10 14 


110/120 


2.700 


0.641 


BH 


BH 


BH 


BH 


F3 


7.5 x 10 14 


100/120 


2.723 


0.858 


BH 


BH 


BH 


BH 


F4 


7.5 x 10 14 


90/120 


2.738 


0.972 


BH 


BH 


no BH 


BH 


F5 


7.5 x 10 14 


80/120 


2.742 


1.010 


BH 


BH 


no BH 


BH 


G4 


7.75 x 10 i4 


90/120 


2.874 


0.957 


BH 


BH 


BH 




G5 


7.75 x 10 14 


80/120 


2.879 


0.993 


BH 


BH 


BH 




HO 


8.0 x 10 14 


1 


2.940 


0.0 


BH 


BH 


BH 


BH 


H5 


8.0 x 10 14 


80/120 


3.016 


0.978 


BH 


BH 


BH 


BH 



Adopted equations of state 
b c 



TABLE IV. Summary of the outcome in the iron core collapse with rigidly rotating initial models. The adopted equations of 
state are discriminated by a single alphabet: 'a' for (ri,r 2 ,pi 4 ) = (1.32,2.25,2.0), 'b' for (1.3,2.50,2.0), V for (1.3,2.22, 1.0), 
and 'd' for (1.28, 2.75, 2.0). "BH" implies that a black hole is formed in a dynamical time scale of the collapse, "no BH" implies 
that a black hole is not formed promptly. " — " in the last column implies that we did not perform simulations for such models. 



model 


Ko (cgs) 


A 


Madm/Mq 


q 


outcome 


D20dlb 


7.0 x 10 14 


1.0 


2.437 


0.647 


BH 


D22dlb 


7.0 x 10 14 


1.0 


2.441 


0.698 


BH 


D23dlb 


7.0 x 10 14 


1.0 


2.445 


0.746 


no BH 


D25dlb 


7.0 x 10 14 


1.0 


2.450 


0.791 


no BH 


D15d05b 


7.0 x 10 14 


0.5 


2.424 


0.440 


BH 


D17d05b 


7.0 x 10 14 


0.5 


2.430 


0.541 


BH 


D20d05b 


7.0 x 10 14 


0.5 


2.437 


0.626 


no BH 


D25d05b 


7.0 x 10 14 


0.5 


2.444 


0.702 


no BH 


H5dl 


8.0 x 10 14 


1.0 


3.019 


0.978 


BH 


H5d05 


8.0 x 10 14 


0.5 


3.020 


0.941 


no BH 



TABLE V. Summary of the final outcome in the iron core collapse with differentially rotating initial models. The adopted 
equation of state is type 'b' for all the cases; i.e., (Fi, T2, P14) = (1-3, 2.5, 2.0). 
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FIG. 6. The distribution map of the outcome in the iron core collapse for rigidly rotating initial models. The horizontal and 
vertical axes denote the spin parameter q and Madm- The filled squares denote the models whose final outcome is a black hole 
irrespective of the equations of state. The open stars denote the models whose final outcome is a black hole for the equations 
of state of 'a' and 'b', and a neutron star for the equation of state 'c'. The crosses denote the models whose final outcome is a 
black hole only for the equations of state 'a'. The open squares indicate no black hole formation irrespective of the equations 
of state. The solid, dotted, and dashed curves indicate the threshold mass above which a black hole is formed as a result of 
the collapse for the equations of state 'a', 'b', and 'c', respectively (see Sec. IVB4 for locating these curves). The long dashed 
curve denotes the mass-shedding limit for rigidly rotating initial models; i.e., no equilibrium configurations for rigidly rotating 
initial condition exist in the right hand side of this curve. 



hole is not formed promptly (models "no BH"), a longterm fallback of matter may lead to black hole formation. In 
this paper, we do not study such longterm black hole formation. Since the simulations are performed for the iron core 
(only part of the whole star), we stop the simulations when the shock front reaches its surface located at a radius 
of <~ 1000 km. To follow the collapse due to the fallback, it would be necessary to take into account not only the 
iron core but the carbon-oxygen envelope. Also, neutrino cooling that is ignored in this paper will play a role for a 
longterm fallback with the duration > 100 ms [71]. 

Figure 6 illustrates that the threshold mass for the prompt black hole formation depends on the adopted equation 
of state as well as the angular momentum. In the following section, we describe, in more detail, the dependence of 
the threshold mass on thermal pressure, rotation, differential rotation, and adopted equations of state, separately. 



1. Contribution of thermal pressure 



First, we focus on spherical collapse. Because of no rotational effect in this case, the threshold mass depends only 
on equations of state. Figure 6 shows that the threshold mass for the prompt black hole formation locates between 
~ 2.1M Q (model AO) and s=s 2.3M Q (model CO) depending on the adopted equations of state. This value is by 20-40% 
larger than the maximum allowed mass for the cold spherical polytrope « 1.6M Q . 

For the equations of state 'a' and 'b', the threshold mass is w 2.1-2.2M Q . On the other hand, for the equation of 
state V in which the value of Y\ is the same as the equation of state 'b' but the value of p nuc is smaller, the threshold 
mass is « 2.3M Q . This fact suggests that magnitude of the thermal pressure P t h generated by the shocks is larger 
for the equations of state with the smaller value of p nuc . This difference in the strength of the shock results from the 
difference in the collapse dynamics. Detailed discussions about the dynamics of the collapse are presented in the next 
subsection. 
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FIG. 7. Snapshots of the density contour curves and velocity vectors in the £-z plane for model D2 with the equation of 
state 'a', at t = 193.29, 193.58, and 193.77 ms. The density contour curves are drawn for p/p max = lCT ^ , (j = 0, 1, 2, • • • 20) 
where p m ax is the maximum density at the selected time slices. The thick solid curve in the last panel denotes the location of 
the apparent horizon. 
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FIG. 8. The same as Fig. 7 but for model D2 with the equation of state 'b' at t = 113.12, 116.38, and 117.10 ms. The thick 
solid curve in the last panel denotes the location of the apparent horizon. 



2. Dependence on Ti and T2 

As discussed and illustrated in [29,34,36], the most important parameter for the dynamics of the collapse during 
the infall phase is Ti. For the smaller value of Ti (for the larger value of |Ti — 4/3|), the depleted pressure at t = is 
larger. As a result, the collapse is accelerated more and the elapsed time in the infall phase is shorter. Furthermore, 
since the depleted fraction of the pressure is larger in the central region than in the outer region for the smaller value 
of Ti, the collapse in the central region proceeds more rapidly. Accordingly, the iron core contracts less coherently 
and only the collapse in the inner region of the iron core is accelerated. Therefore, the mass of the inner core at the 
bounce phase is smaller for the smaller value of IV The smaller mass of inner core indicates that the fraction of the 
iron core which undergoes the shock heating is larger since the shock wave is generated at the outer edge of the inner 
core. Therefore, the contribution of the thermal pressure to the inner core is more important for models with the 
smaller value of IV 

Since the mass of the inner core in the infall and bounce phases is larger for the larger value of IV the degree of 
the overshooting of the inner core at the bounce phase is larger, due to the stronger gravitational attraction force. 
Also, for the smaller value of Y2, the degree of the overshooting at the bounce phase is larger since a relatively higher 
density is required to supply the sufficient pressure to halt the collapse. In the case that mass of the inner core is 
not so large as to collapse to a black hole, the larger overshooting results in stronger shock generation, since the 
stored energy of the inner core at the bounce becomes larger (cf. Sec. IV D 1). However, the larger overshooting also 
increases the risk of collapsing to a black hole since the compactness (the value of a c ) at the bounce becomes larger 
(smaller). 

To understand the dependence of the collapse on equations of state in more details, we generate Figs. 7 and 8, 
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FIG. 9. The same as Fig. 7 but for model F5 with the equation of state 'a', at t = 198.74, 199.03, and 199.35 ms. The thick 
solid curve in the last panel denotes the location of the apparent horizon. 



in which we display the snapshots of the density contour curves and velocity vectors in the x-z plane for models 
D2a and D2b. For models with the equation of state 'a' in which Ti = 1.32 and T2 — 2.25, the collapse proceeds 
coherently due to the fact that T\ is close to 4/3, and hence, the mass of the inner core increases to be larger than 
the maximum allowed mass of the corresponding rotating neutron star at the bounce. Consequently, the inner core 
promptly collapses to a black hole without generating shocks (cf. Fig. 7). This result indicates that for the smaller 
value of |Ti — 4/3 1 , the prompt black hole formation is more liable and the threshold mass of the iron core for the 
prompt black hole formation is smaller. 

For the equations of state 'b' and 'd', on the other hand, the collapse proceeds less coherently due to the smaller 
values of Ti (the larger values of |Ti — 4/3 1 ) . As a result, the mass of the inner core at the bounce is below the 
maximum allowed mass, and hence, the inner core experiences a bounce and forms a compact protoneutron star for 
a while before collapsing to a black hole. Also, the shock is formed at the outer edge of the inner core and then 
propagates outward (cf. Fig. 8). This implies that the thermal pressure which helps supporting the inner core is 
generated in contrast to the case for model D2a. Due to this difference, the threshold mass of the iron core for the 
prompt black hole formation with the equation of state 'b' is larger than that for 'a'. However, the shock formed is 
not strong enough to blow up the sufficient matter outward even for model D2b: A part of the matter in the outer 
region falls into the inner core located at the center. This fall-back leads to an accretion induced collapse to a black 
hole (see the third panel Fig. 8). 

On the other hand, we do not find significant difference in the threshold mass for the prompt black hole formation 
between the equations of state 'b' and 'd' (see table IV). The plausible reason is described as follows. First, since the 
mass of the inner core formed at the bounce is smaller for models with the equation of state 'd', the larger fraction of 
the matter inside the protoneutron star will experience the shock heating. On the other hand, the shock wave itself 
is weaker for the equation of state 'd' since the amplitude of the core bounce and the value of T t h (= Ti) is smaller. 
These two contrary effects cancel each other resulting in approximately the same threshold mass for the prompt black 
hole formation. A more sophisticated parameter study is required to clarify the small difference in the threshold mass. 



3. Dependence on p nuc 



In Fig. 4, we show the evolution of p c and a c for models Die and D15c, for which a black hole is formed promptly 
as a result of the collapse. For comparison, the result for model Dlb is shown together. The evolution of p c and a c 
with two equations of state 'b' and 'c' are identical during the infall phase. However, in the bounce and subsequent 
phases, increase of p c (decrease of a c ) is delayed and the value of p c (a c ) is smaller (larger) for the equation of state 
'c'. This is a result from the fact that the sudden stiffening occurs at an earlier stage of the collapse with the equation 
of state 'c' due to the smaller value of p nuc - 

Figure 4 also indicates that the amplitude of oscillation of the inner core at the bounce is larger for models with 
the equation of state 'c'. Several effects are responsible for this result. First, the small value of p nuc in the equation of 
state V results in a smaller mass of the inner core and smaller averaged density at the bounce. As mentioned in Sec. 
IV B 2, the smaller mass of the inner core results in the fact that a larger amount of the matter experiences the shock 
heating. The smaller density implies that the pressure at the surface of the inner core is smaller, and hence, the work 
exhausted by the inner core in converting the oscillation energy to the infalling outer envelop is smaller. Accordingly, 
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FIG. 10. The same as Fig. 7 but for model F5 with the equation of state 'b', at t = 114.50, 116.38, 119.19, 124.90, and 
131.86 ms. The thick solid curve in the last panel denotes the location of the apparent horizon. 



the amplitudes of the bounce and the subsequent oscillation of the inner core become larger. Also, the equation of 
state V is "stiffer" than 'b' in the (higher) density range between p w 1.0 x 10 14 and « 2 x 10 15 g/cm 3 (cf. Fig. 1). 
This also contributes to generating stronger shock waves. Due to these effects, the smaller value of p nuc results in the 
larger threshold mass for the prompt black hole formation. 



4- Effects of the rotation 

As the value of q increases, the threshold mass for the prompt black hole formation becomes larger because the effect 
of the rotation effectively supplies additional pressure to the iron core and reduces the amount of matter falling into 
the central region. The threshold mass may be written approximately as a quadratic form [44]: M t hreshoid = k q 2 + Mo, 
where Mo is the threshold mass for the spherical case. This reflects that the rotational kinetic energy (or the centrifugal 
force) depends on q 2 . Assuming this form of the threshold mass, we approximately draw the threshold curves in Fig. 
6. It is found that the rotational effect increases the threshold mass by < 15-25% for q ~ 1. This value is comparable 
with the amplification factor of the maximum mass for rigidly rotating neutron stars in equilibrium [67] . 

The coefficient, k, depends on the strength of the shock and the adopted equation of state. The threshold mass for 
the prompt black hole formation at the maximum value of q is larger for the models with the smaller value of p nuc 
as discussed in the previous subsection (compare models with equations of state 'b' and 'c' in Table IV). However, 
the dependence of the threshold mass at the maximum value of q on Ti is not very simple. Comparing models with 
equations of state 'a' and 'b', the black hole is more liable to be formed with equation of state 'a', i.e., for the larger 
value of Ti. On the other hand, comparing models with equations of state 'b' and 'd', no difference is found in 
the present numerical study. This is because the process of the black hole formation for models with Ti = 1.32 is 
significantly different from those with Ti = 1.30 and 1.28, as described in Sec. IV B 2. 

To sec the effect of the rotation on the collapse dynamics in more detail, in Figs. 9 and 10, we display the snapshots 
of the density contour curves and velocity vectors at selected time slices for models F5a and F5b. For these models, 
the value of q is nearly maximum among the rigidly rotating initial models of a given mass. For model F5a, the 
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(a) (b) 

FIG. 11. The spin parameter distribution inside differentially rotating iron cores of (a) D25dl (solid curve) and D22dl 
(dashed curve); (b) D20d05 (solid curve) and D17d05 (dashed curve). The two dotted curves in both figures denote the spin 
parameter distribution for the rigidly rotating models D3 (upper dotted curve) and D25 (lower dotted curve). Note that for 
model D3b, a black hole is formed promptly, while for D25b, a neutron star is formed. 



collapse proceeds more rapidly in the direction of rotational axis (z axis) than in the equatorial plane during the infall 
phase. This is because the centrifugal force is stronger in the equatorial plane than around the z axis (see the first 
panel of Fig. 9). Accordingly, the collapsing inner region is deformed to be an oblate shape. As the collapse proceeds, 
the density contour in the outer part of the inner core is deformed to be a concave structure. At the bounce, this leads 
to formation of a steep density gradient along the rotational axis (see the second panel of Fig. 9). This would result 
in stronger shock waves along the rotational axis than in the equatorial plane. For model F5a, however, the mass of 
the inner core at the bounce is so large that the strong gravitational attraction force prevents the shock from being 
propagated. Consequently, the inner core promptly collapses to a black hole. On the other hand, for less massive 
cases, shock waves propagate outward. In the shock propagation, this asymmetry generates anisotropic shocks (see 
also Sec. IV D). 

For model F5b, shock waves propagate outward since the mass of the inner core at the bounce is smaller than the 
maximum allowed mass due to the smaller value of T\. In this case, however, the density gradient along the z axis 
is not so steep as that for model F5a (see the first panel of Fig. 10), although the mass and the angular momentum 
of the progenitor arc identical between F5a and F5b. As a result, the asymmetry in the shock front is small (see the 
second to fourth panels of figure 10). The reason is that the mass of the inner core at the bounce is smaller for model 
F5b. For rigidly rotating initial models, the centrifugal force is stronger for larger cylindrical radius, and hence, the 
smaller mass of the inner core implies that the effect of rotation is less important. 

For model F5b, the shock is strong enough to reach the surface of the iron core (see the fifth panel of Fig. 10). 
However, in a region behind the shock, the fluid elements (in particular) around the rotational axis fall back into a 
protoneutron star formed at the center, and eventually, the protoneutron star collapses to a black hole (see the last 
panel in Fig. 10). This illustrates that a black hole may be formed in a rather long time scale even for models in 
which the black hole is not formed promptly, if the progenitor is rotating. 



5. Effects of differential rotation 

Table V shows the outcomes of the collapse for all of differentially rotating initial models listed in Table III. In 
the simulation for differentially rotating models, we only adopt the equation of state 'b'. Model D25dl with A = 1.0 
does not form a black hole, while its rigidly rotating counterpart D25 for which the mass and the angular momentum 
are approximately the same as those of D25dl collapses to a black hole. For differentially rotating model D20d05 for 
which A = 0.5, and the mass and the angular momentum arc slightly smaller than those for rigidly rotating model 
D2, black hole is not formed. For A = 0.5, black hole is not formed promptly even from a very massive initial model 
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H in which Madm ~ 3.0M Q (cf. Table V). These indicate that the black hole is less liable to be formed promptly for 
the higher degree of differential rotation. 

The results found here are quite reasonable since with the decrease of A, the angular velocity in the inner region 
increases. To clarify this effect in a more quantitative manner, we generate Fig. 11, which shows the spin parameter 
distribution q*{j) defined by Eq. (14). First, note that a black hole is formed for models D25 and D22dl but not 
for D3 and D25dl with the equation of state 'b'. The two dotted curves in Figs. 11(a) and 11(b) denote q*(j) 
of these two models D3b (upper dotted curve) and D25b (lower dotted curve), respectively. The solid and dashed 
curves in Fig. 11(a) denote q*(j) for models D25dl and D22dl for which A = 1.0. Although the total value of the 
spin parameter of the differentially rotating iron cores of D25dl (q — 0.791) is smaller than that for rigidly rotating 
model D3 (q — 0.888), the value of q* around the central region (m*(j) < M*/2) exceeds that of D3. As shown in 
[12] in detail, the value of q* around the center plays a crucial role in determining the criterion for prompt black 
hole formation; i.e, the iron core with a larger value of g» around the center is less liable to collapse to a black hole. 
Therefore, it is reasonable that black hole is not formed for model D25dl. On the other hand, the value of around 
the center for model D22dl (q = 0.698) is as large as that for model D25 around the central region. This implies the 
rotational centrifugal force is not strong enough to prevent the iron core from collapsing to a black hole. 

Figure 11(b) shows q*(j) for models D20d05 (solid curve) and D17d05 (dashed curve) with A — 0.5. A black hole 
is not formed for D20d025 while it is formed for D17d05. Figure 11(b) clarifies that the spin parameter distribution 
is flatter with the smaller value of A. Therefore, the values of the global spin parameter q for A = 0.5 are much 
smaller than those of the rigidly rotating counterparts. On the other hand, the value of q* around the central region 
becomes larger. As a result, black hole formation is more effectively prevented. This quantitatively indicates that the 
threshold mass for the prompt black hole formation is larger for models with larger degree of differential rotation. 



C. Prediction of the final system 




0.6 0.7 0.8 0.9 1 

m,(j)/M + 

FIG. 12. The distribution of q(j) (top panel) and the value of jisco as a function of m*(j)/M* for models H5 (solid curve), 
F5 (dashed curve), and D3 (dotted curve). The dotted dashed vertical line denotes m»(j)/M* = 0.9 where jisco for F5 takes 
the maximum value approximately. 

Because of the assumption that the viscous effect is negligible during the collapse, the specific angular momentum 
j = it v of each fluid clement is conserved in the axisymmctric system and so do the functions m„{j) and J*(j). Using 
this fact, the final outcome after the black hole formation can be predicted [19,10,12]. 

Let us consider innermost stable circular orbit (ISCO) around the growing black hole located at the center. Assuming 
that a fluid element around a black hole is approximately in a circular orbit, it will fall into the seed black hole 
eventually if its value of j is smaller than that at the ISCO (jisco)- The value of jisco depends on the mass and 
the angular momentum of the black hole and changes as the ambient fluid elements accrete onto the black hole. If 
Jisco increases as a result of the accretion, the more ambient fluid elements will fall into the black hole. On the 
other hand, if jisco decreases during the accretion, no more fluid element will fall into the black hole, and as a result, 
the dynamical growth of the black hole will terminate. Therefore, if jisco has a maximum (hereafter denoted as 
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JiSCOrmax), the black hole will grow until j reaches jiscOrmax- Namely, it is reasonable to expect that the final values 
of the mass and the spin parameter of the black hole will be 

M B H « ™-*(jlSCO:ma X ), (36) 
<7BH ~ <7*(jlSCO:ma X )- (37) 

with the disk mass M A dm - ™*(jiscO:max)- 

To estimate the value of jisco , we assume that the spacetime metric can be instantaneously approximated by that 
of a Kerr spacetime of the mass m*(j) and the spin g*(j). On these approximations, we can compute jisco 01 a 
growing black hole as [72,13], 



\/m*(j)nsco (nsco 2 - 2q*(j) m *ti)V m *ti) r isco + (g* (j')ttt.* (j)) 2 ) 

JISCO = " — 1 ( 38 ) 

nsco (nsco 2 - 3m*(j)nsco + 2g* (j)m* m* (j)nsco 

where 

nsco = m*(i) [3 + ^2 - {(3 - Zi)(3 + Zi + 2Z 2 )} 1/2 
Z 1 = l+[1- cUj) 2 ] 1/3 [{1 + g*(i)} 1/3 + {1 - 9 ,(i)} 1/3 
Z 2 =[3g.(j) 2 + Z 1 2 ] 1/2 . 

In Fig. 12, we show the quasi-local spin parameter distributions inside the iron core, defined by Eq. (14), and the 
value of jisco evaluated by Eq. (38) for models H5, F5, and D3. As described above, the dynamical evolution of the 
formed black hole will terminate when jisco reaches a local maximum. Then, the mass, Mbh, and spin, qbh, of the 
the dynamically evolved black hole are given by Eqs. (36) and (37). 

Figure 12 indicates that for model F5 and H5, the inner region collapse to form a black hole of Mbh/M* « 0.90 
and qbh ~ 0.76. On the other hand, the matter in the outer region of high specific angular momentum will form a 
massive disk of M^k/M* w 0.10 around the black hole. Similarly, for model D3, a black hole of Mbh/M* w 0.93 and 
<7bh ~ 0.73 will be formed. Such a rapidly rotating black hole surrounded by a massive disk is one of the promising 
candidates for the central engine of the long-duration gamma-ray bursts [7-9] . 

Although the above prediction is quite reasonable [12], confirmation of this prediction requires to carry out a 
simulation until the black hole plus disk system is formed. To accomplish this, the so-called black hole excision 
techniques are required. For the case of the collapse with stiff (r = 2.0) equations of state, we have confirmed it 
possible to continue the simulation more than ~ 100M after the formation of black hole by using techniques based on 
the so-called simple excision developed by Alcubierre and Brugmann [73] . We plan to wrestle with the confirmation 
of the above prediction by extending the techniques to the collapse with soft equations of state. 



D. Protoneutron star formation 



1. Dependence of explosion on the equations of state 



As discussed in the previous subsections, the dynamics of the collapse and the threshold mass for the prompt black 
hole formation depend sensitively on the equations of state. In the case of neutron star formation, we also find that 
features of the explosion depend sensitively on the equations of state. We illustrate this fact focusing on rapidly 
rotating model D5 in which black hole is not formed promptly irrespective of the equations of state. 

In Fig. 13, we display the contour plots for model D5a. For the collapse with the equation of state 'a' which is 
stiff in the subnuclear density (with larger value of Y\) while soft in the supranuclear density (with smaller value of 
^2), the steep density gradient is formed around the rotational axis of the inner core. On the other hand, such steep 
gradient is not formed around the equatorial plane (see the second panel of Fig. 13). Due to the steep gradient, strong 
shock waves are generated along the rotational axis. Consequently, a jet-like explosion is seen (see the third panel 
of Fig. 13). The explosion velocity in the direction of the rotational axis is much larger than that in the equatorial 
plane (see the fourth and fifth panels of Fig. 13): The maximum speed becomes « 0.8c around the rotational axis 
near the shock. Thus, the shock wave reaches the surface of the iron core much more quickly (see the last panel of 
Fig. 13). At this time, a funnel structure is formed around the rotational axis. 
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FIG. 13. The same as Fig. 7 but for model D5 with the equation of state 'a', at t = 199.47, 199.69, 200.31, 202.01, 203.33, 
and 205.84 ms. 
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FIG. 14. The same as Fig. 7 but for model A5 with the equation of state 'a', at t = 201.53, 201.88, and 203.69 ms. 
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FIG. 15. The same figure as Fig. 7 but for model D5 with the equation of state 'b', at t = 114.72, 115.63, 116.75, 118.60, 
120.37, and 124.25 ms. 



The structure of bipolar shock waves depends on mass of progenitor. For model A5a in which Madm ~ 1.97Mq, a 
steep density gradient is also formed along the rotational axis (see Fig. 14), but it is formed only in narrower region 
than that for model D5a (compare the first panel of Fig. 14 and the second panel of Fig. 13). Reflecting this result, 
the structure of the shock front is sharper than that for model D5a and is deformed to be spear-like (compare the 
second panel of Fig. 14 and the third panel of Fig. 13). 

Such bipolar and energetic supernova explosion is required to explain the observations of several hypernovae, which 
are suggested to be associated with a gamma-ray burst [6,74-77]. Indeed, the energetic explosion along the rotational 
axis is preferable to avoid the baryon- loading problem in the fire-ball model of the gamma-ray bursts [9] . The present 
results suggest that the collapse of rapidly rotating iron cores of M > 2Mq with a class of equations of state similar 
to 'a' is preferable for realizing such a state. 

To see the dependence of the features of the explosion on the equations of state, we also display the contour plots 
for models D5b and D5d in Figs. 15 and 16. As the first panel of these figures show, the magnitude of the density 
gradient along the rotational axis is not very different from that in the equatorial plane, and hence, structure of the 
shock is not very aspherical (see the second panel of the Figs. 15 and 16). 

The shock for model D5b is strong enough to quickly propagate through the outer core, sweeping the matter falling 
into the shock. For model D5d, on the other hand, the shock is not as strong as that for D5b, and hence, the 
propagation speed is slow. Also, the density gradient formed around the shock is not very different between along 
the rotational axis and in the equatorial plane, and hence, the structure of the shock becomes eventually oblate (not 
prolate) due to the effect of the rotation (see the third panel of Fig. 16). As the matter behind the shock falls and 
beats the protoneutron star at the center, it oscillates and the shocks are formed continuously. Thus, the shocks 
propagates outward and helps sweeping the infalling matter outward (see the fourth and fifth panels of Fig. 16). The 
features found for models D5b and D5d are not found in model D5a. This fact implies that the explosion mechanism 
depends sensitively on the equations of state for rapidly rotating iron core collapse. 

It should be also mentioned that the bipolar explosion has not been found in the previous study [36] in which the 
same parametric equation of state with (Ti, T 2 ) = (1.32, 2.5) and smaller core mass of 1.5M Q is adopted. (Note that 
this equation of state is stiffer in supranuclear density than the equation of state 'a' of the present paper). This fact 



22 



t=63 93 (ms) 

I I I I I I I / / / / / /- 

l~T~l~r~t~i—LJJ / / Pi 

I i I I I I J />-< / / 




60 



4 



t=84 76 (ms) 



I I I J / / / / 
//////// 
IJi//////// 

liy//////// 
i j i i / / 

1 2 ////// ✓ 
1 1 1 J J / s s s 
2 / JJ/sssss 





60 



t.-B6 78 (ms) 



■12222////// p?^ 

1 i 2 /// 
J J J////// 




t- 90 24 (ms) 

TTTTTTryTyyy 
1 1 1, 1 1 1 j j / / / £ 



TTTT7777V-JS . 
J]////////// 
' '///// 
' / / / 



,3? 




20 40 60 
X (km) 



80 



t= 90 34 (ms) 

TTTTTTTyyyT 
I I I 1 1 1 J J / / / »£ 



I I I J J JTTV~^~> _ 
12 2 2//////// 



/ ? iz 




40 60 

X (km) 



100 



50 



1=100.01 (ms) | 

I i / j' i '/'/)///'/ '/ 77>S 
iiJi//////// / z^/- 

kiii/y////// / «o«3<?- 




50 100 
X (km) 



FIG. 16. The same as Fig. 7 but for model D5 with the equation of state 'd\ at i = 83.92, 84.76, 86.78, 88.92, 90.24, and 
100.01 ms. 



implies that an equation of state with not only a large value of Fi(= 1.32) but also a small value of T 2 < 2.25 may 
be required for producing the bipolar explosion. 



2. Rotational profile of the protoneutron star 

In this subsection, we present the rotational profiles of the inner region during the collapse and consider the 
possibility for the onset of the nonaxisymmetric instabilities. Figures 17-20 show the angular velocity profile along 
the x axis at selected time slices for models D5a, D5b, D5c, and D5d, respectively. The first panels of Figs. 17-20 
show the angular velocity profile around the bounce phase. As the collapse proceeds, the angular velocity around 
the central region increases toward the maximum value, f2 max , achieved at the bounce. The values are f2 max ~ 8000, 
10000, 6000, and 7500 s _1 for models D5a, D5b, D5c, and D5d, respectively (see dotted lines in the first panels). 
The inner core at the bounce, when the maximum angular velocity is achieved, is approximately rigidly rotating 
irrespective of the equations of state, although the fluid elements outside the inner core is differentially rotating as 
a whole. These results indicate that the differential rotation is not enhanced in the inner core during the infall and 
bounce phases. 

The second panels of Figs. 17-20 show the angular velocity profiles in the ringdown phase, from which we find that 
differential rotation is enhanced in the ringdown phase due to the oscillation of the protoneutron star and the infall of 
the fluid elements of high specific angular momentum from the outer region. However, the degree of the differential 
rotation is not very high. This feature is also found in a Newtonian simulation [33]. These panels also indicate that 
the amplitude of the oscillation of the angular velocity is larger for models with the equation of state 'c'. This is due 
to the fact that the amplitude of the oscillation of the protoneutron star is larger as described in Sec. IV B 3. 

In the post-bounce oscillation phase, the oscillation of the protoneutron star gradually damps. The damping is 
caused by the compressional work of the protoneutron star to the infalling matter from the outer envelop. As a result 
of damping, a quasistationary protoneutron star rotating rapidly at the period P = 2n/Q <~ 1 ms is eventually formed. 
Note that inner region of the formed protoneutron star is rigidly rotating (see the third panel of Fig. 18). 
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FIG. 17. The angular velocity profile inside the iron core along the :r-axis for model D5a at selected time slices: (a) at the 
bounce phase of t = 199.43 (solid curve), 199.62 (long dashed curve), 199.73 (dashed curve), and 199.86 ms (dotted curve); (b) 
at the ringdown phase of t = 206.49 (solid curve), 206.75 (long dashed curve), 206.94 (dashed curve), and 207.12 ms (dotted 
curve) . 




FIG. 18. The same as Fig. 17 but for model D5b at selected time slices: (a) at the bounce phase of t = 114.39 (solid curve), 
114.61 (long dashed curve), 114.74 (dashed curve), and 114.87 ms (dotted curve); (b) at the early ringdown phase of t = 118.86 
(solid curve), 119.16 (long dashed curve), 119.29 (dashed curve), and 119.53 ms (dotted curve); (c) at the late ringdown phase 
of t = 122.06 (solid curve), 122.36 (long dashed curve), 122.67 ms (dashed curve). 




FIG. 19. The same as Fig. 17 but for model D5c at selected time slices: (a) at the bounce phase of t = 114.46 (solid curve), 
114.59 (long dashed curve), 114.72 (dashed curve), and 114.83 ms (dotted curve); (b) at the ringdown phase of t = 120.96 
(solid curve), 121.11 (long dashed curve), 121.67 (dashed curve), and 121.97 ms (dotted curve). 
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FIG. 20. The same as Fig. 17 but for model D5d at selected time slices: (a) at the bounce phase of t = 83.78 (solid curve), 
83.91 (long dashed curve), 84.04 (dashed curve), and 84.15 ms (dotted curve); (b) at the ringdown phase of t = 86.10 (solid 
curve), 86.32 (long dashed curve), 86.51 (dashed curve), and 86.71 ms (dotted curve). 



To infer the nonaxisymmetric stabilities, we estimate the value of T rot /W at the bounce. Note that in general 
relativity, there is no unique definition for T TOt /W for dynamical spacetimes. In this paper, we define the rotational 
kinetic energy T rot , and the gravitational potential energy W, for dynamical spacetime by 

T rot = i J d 3 xp,u v v^, (39) 

W = Ttotai + U, where U = J d 3 xp*e. (40) 

To define W, we use the fact that Madm ~ M* — W + Ttotai + U and M* « Madm- Total denotes a total kinetic 
energy. Unfortunately, we do not know the appropriate definition for Total- However, as far as configurations at the 
maximum compression and at a final relaxed state are concerned, the rotational kinetic energy T rot is nearly equal to 
Total, an d thus, Ttotai m Eq. (40) may be replaced by T rot . In this definition, T TOt /W will give slightly overestimated 
values in other stages such as the infall phase. The value of T TOt /W at the bounce is « 0.10 and « 0.13 for models D5b 
and D5a. These values are much smaller than the plausible critical value of T rot /W for the onset of the dynamical 
instability T TOt /W ~ 0.24-0.25 for nearly rigidly rotating stars [68,78,65], and hence, the inner core at the bounce 
and formed protoneutron star is unlikely to be unstable against nonaxisymmetric perturbations in the dynamical 
time scale. However, these values are as large as the critical value of T m t/W for the secular instabilities driven by 
gravitational wave emission [79,80,65], or driven by viscosity [81,65]. Therefore, the formed protoneutron star may be 
unstable against nonaxisymmetric perturbations in secular time scales 3> 100 ms. 

A sufficiently rapidly and differentially rotating iron core may collapse to form a protoneutron star of T TOt /W > 0.25 
and become unstable against nonaxisymmetric perturbations. Such a rotating nonaxisymmetric object will be a strong 
gravitational wave emitter. In a companion paper [38] , we study conditions for the onset of the dynamical instabilities 
in three-dimensional numerical simulation in full general relativity. We find that the dynamical instabilities set in 
only for the case that the progenitor is highly differentially rotating with A < 0.1. Thus, the assumption of axial 
symmetry in this paper is justified. 



E. Gravitational waves 

1. Gravitational waveforms 

Gravitational waveforms are computed in terms of the quadrupole formula described in Sec. II C. As illustrated 
in a previous paper [55], approximate gravitational waveforms can be computed even for highly relativistic, highly 
oscillating, and rapidly rotating neutron stars using the quadrupole formula except for systematic underestimate of 
the amplitudes of 0{M/R). In the case of protoneutron star formation, gravitational waves are emitted mainly by 
its oscillation. Thus, with the quadrupole formula, it is possible to approximately compute gravitational waveforms 
emitted during rotating iron core collapse to a protoneutron star. 

In Figs. 21-23, gravitational waves for models D2-D5 with equations of state 'a', 'b', and 'd' are shown. The 
gravitational waveforms during the black hole formation evaluated by the quadrupole formula are also presented 
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FIG. 21. Gravitational waveforms computed by the quadrupole formula described in Sec. II C for models D5a (solid curve) 
and D4a (dotted dashed curve), and D3a (dashed curve). 



together for models D3a, D2b, and D2d. We note that from A 2 = I zz — I xx , the amplitude of gravitational waves at 
a distance of r from the source is calculated by 

* « 3 >< ltr2 ° ( 77^~) ( — 1 sin2 °> ( 41 ) 
V 1000 cm J V r J ' 

where 8 denotes the angle between the line of sight and the rotational axis, and (sin 2 6>) = 2/3. Thus, the typical 
amplitude is <~ 2 x 10~ 20 for a event at the galactic center. 

Figures 21-23 show that the strong dependence of the dynamics of the collapse on the values of Ti is reflected in the 
amplitude of gravitational waveforms. The following features are in particular worthy to note. First, the amplitudes 
of gravitational waves during the infall and bounce phases are smaller for the smaller value of Ti (the larger value of 
Ti — 4/3 1 ) . The reason is that the smaller value of Ti results in the smaller mass of the inner core at the bounce (cf. 
discussion in Sec. IV B 2). This feature agrees with that found in previous works [29,34,36]. 

Second, the time-averaged amplitudes of the gravitational waves in the ringdown phase are negative for models with 
equations of state 'b' and 'd' while positive for the model with 'a'. This results from the difference of the dynamics 
in the infall phase. For the equations of state 'a', the matter collapses and bounces in a coherent manner, while for 
the equations of state 'b' and 'd', in a less coherent manner. In the case 'b' and 'd', therefore, the matter in the outer 
region infalling toward the center suppresses the oscillation of the protoneutron star. On the other hand, for models 
with 'a', such suppression dose not work effectively because of the smaller mass fraction in the outer region. 

Third, the gravitational waveforms for rapidly rotating models with the equations of state 'a' (models D4a and D5a 
in Fig. 21) are qualitatively different from those for others. The modulation of the gravitational waveforms are quite 
remarkable. This reflects the bipolar explosion, in which absolute value of I zz is much larger than that of I xx . 

Fourth, the so-called "type-Ill" gravitational waveforms [29] are not found in any model. Previous studies have 
indicated that type-Ill waveforms are generated for a small value of Y\ ~ 1.28 [29,34,36]. In these waveforms, the 
amplitude of the first spike is significantly suppressed to the value \I ZZ — I xx \ ~ 30 cm. This is because the mass of 
the protoneutron star at the bounce phase is very small for such small value of T\. In the present case, the mass 
of the adopted iron core is much larger than that for the previous studies [29,34,36]. As a result, the mass of the 
protoneutron star is not very small, and hence, the type-Ill gravitational waveforms are not generated. 

For a given mass, the amplitude of gravitational waves is increased with the increase of the angular momentum in 
the present numerical results. However, this is not trivial result because of the following reasons. First, recall that the 
amplitude of gravitational waves is proportional to the quadrupole moment (\I ZZ — I xx \) and the inverse square of the 
dynamical time scales, T c h a r, of the system. The value of \I ZZ —I xx \ increases as the total angular momentum increases 
since the radius and the degree of deformation of the inner core become larger due to the increasing centrifugal force. 
On the other hand, the characteristic time scale becomes longer as the angular momentum of the inner core increases 
since the rotation effectively supplies additional pressure to the inner core and the collapse is decelerated. Thus, these 
two contrary effects may cancel each other with a certain value of spin parameter q, resulting in the saturation of 
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FIG. 22. The same as Fig. 21 but for models D5b (solid curve), D3b (dotted dashed curve), and D2b (dashed curve). 




FIG. 23. The same as Fig. 21 but for models D5d (solid curve), D3d (dotted dashed curve), and D2d (dashed curve). 



the burst-like gravitational wave signals. Indeed, such saturation was found at q ~ 1 in previous Newtonian studies 
(e.g., [28]) since at such a large value of q, the iron core does not collapse to a sufficiently compact state. However, 
saturation is not found in the general relativistic studies [34,36]. This implies that the general relativistic gravity is 
strong enough to overcome the centrifugal force for the rigidly rotating case with q < 1. This effect also suppresses 
to yield type-II gravitational waves [34] which appear in the case that the centrifugal force is too strong for the iron 
core to reach the nuclear density. 

In Fig. 24, we show gravitational waveforms for models B5b, D5b, and F5b to illustrate the dependence of them 
on the mass of the iron core. In the infall phase, the amplitude of gravitational waves is larger for the larger-mass 
model. This is quite natural since for the larger-mass model, the quadrupolc moment of the inner cores are larger. 
On the other hand, the amplitude of the strong spike at the bounce saturates when the mass of the iron core reaches 
w 2.5M Q . The plausible reason is that the inner core shrinks to be so compact that the quadrupole moment is 
decreased significantly 

The effect of mass on the amplitude of gravitational waves can be analyzed by comparing the present results with the 
previous results in [36], in which collapse of the iron core of mass M « 1.5M Q and rotational energy T rot /W ~ 0.009 
are studied. In [36], we found that the amplitudes of gravitational waves at the bounce are \I ZZ — I xx \ « 650,200, 
and 30 cm for models with (fi, T 2 ) = (1.32, 2.0), (1.31, 2.5), and (1.28, 2.5). These results should be compared 



27 



1000 



B5b 



1 500 



110 112 



_J I I I I I L_ 

114 116 
T (msec) 



118 120 



FIG. 24. The same as Fig. 21 but for models F5b (solid curve), D5b (dotted dashed curve), and B5b (dashed curve). 




FIG. 25. The same as Fig. 21 but for models D5b (dashed curve) and D5c (solid curve) 



with the present results for D5a, D5b, and D5d, in which M ~ 2.5M© and T TOt /W w 0.009 and the amplitudes of 
gravitational waves are \I ZZ — I xx \ « 1350, 1200, and 750 cm, respectively. These values are larger than those found 
in [36] by factors of w 2, 6, and 25. Thus, with the increase of the mass, the amplitude is increased in a nonlinear 
manner. For Ti = 1.32, the increase factor is small. The reason is that for this case, the inner core at the bounce 
becomes sufficiently large due to the coherence collapse in the infall phase, and the mass is close to the value for the 
saturation as mentioned above. For other cases, the increase factor is much larger. In particular, the increase factor 
for Ti = 1.28 is outstanding. This is associated with the fact that the dynamics in the bounce phase qualitatively 
changes with the increase of mass. For Ti = 1.28 with a small mass M ~ 1.5M Q , the mass of the inner core formed 
at the bounce is very small, and as a result, type-Ill gravitational waves are emitted [34] with the small maximum 
amplitude as \I ZZ — I xx \ ~ 30 cm. On the other hand, for T\ = 1.28 with a large mass M ~ 2.5M Q , the mass of the 
inner core can be sufficiently large due to its strong self-gravity. As a result, type- 1 gravitational waves are emitted 
and the amplitude of gravitational waves is significantly increased. 

An interesting finding in gravitational waveforms for the larger-mass model is that the first spike in the ringdown 
phase emitted just after the strong negative spike at the bounce decreases as the mass of the iron core increases (see 
Fig. 24) . This spike is associated with the outward motion of the inner core and usually positive for the smaller-mass 
models (such as models of Mw 1.5M studied in [36,34,29]). For models of larger mass such as studied in this paper, 
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on the other hand, the post-bounce first spike becomes negative due to suppression of the outward oscillation by the 
infalling matter. The suppression is larger in particular in the direction of the rotational axis because of the absence 
of the centrifugal force. This decreases the oscillation amplitude to be negative. 

Gravitational waveforms for models D5b and D5c are compared to sec the dependence on p n uc in Fig. 25. It 
shows that differences in the dynamics associated with the difference in the value of p nuc are reflected in gravitational 
waveforms. For example, the amplitude at the bounce for D5c is smaller than that for D5b. This is because the 
central density at the bounce is smaller for D5c, and accordingly, the characteristic time scale is longer. On the other 
hand, the amplitude of the post-bounce first spike is larger for D5c. This is because the amplitude of oscillation of 
the inner core is larger for D5c as described in Sec. IV B 3. 

Before closing this subsection, we comment on the frequency of gravitational waves emitted in the case of black hole 
formation. After the formation of a black hole, the ambient matter falls into it. As a result, the so-called quasi-normal 
modes (e.g. [82]) of the black hole will be excited and gravitational waves associated with the spacetime oscillation 
will be emitted. According to Echevcrria [83], scmicmpirical expression for the real part of the quasi- normal mode 
frequency for I = 2 is given by 

M V 1 / F(q) \ 



f « 12 — p^- kHz, (42) 

■' \MqJ V 37/100 y ' V 1 

where 

1-^(1 - q ?' w . (43) 

Equation (42) indicates that the frequency of gravitational waves associated with the quasi-normal mode is ~ 7.5 kHz 
for models D3 and F5 for which the mass and spin parameter of the formed black hole are predicted to be M « 2.3 
2.4M and q w 0.7 (see Sec. IV C). This value is far out of the best-sensitive frequency range for detection of laser 
interferometric gravitational wave detector [84] . 

In reality, the mass of black hole will significantly increase in a longterm evolution since the matter outside the 
iron core will fall into a black hole. After such accretion, the frequency of the quasinormal mode will be decreased 
significantly. Gravitational waves emitted from a high-mass black hole M > 20M Q may have an appropriate frequency 
for detection by the laser interferometric detectors. 



2. Energy power spectra 

In Fig. 26, we show the energy power spectra of gravitational waves with a mode of I = 2 and m = for selected 
models. For all the models, the maximum values of the spectra arc located at / pca k ~ 1-2 kHz, which are consistent 
with the typical bounce time intervals of these models ~ 0.5-1 ms. The value / pca k ~ 1-2 kHz is by a factor of ~ 2 
larger than the previous Newtonian results [29,33]. This reflects the features in the dynamics of the collapse that 
the central density and compactness of the inner core at the bounce in general relativity are larger than those in the 
Newtonian results due to the stronger attraction force [34]. The peak frequency is slightly larger than that found by 
Dimmelmeier et al. [34] . This is due to the fact that the mass of the iron core adopted in this work is larger than 
theirs. For / > 2 kHz, the spectra decline steeply for all the models. This feature is consistent with that in [34,33]. 

The spectrum for model D5a (see Fig. 26(a)) is rather broad for a low frequency region (100 < / < 1000 Hz) and 
flatter than those for other models. Such characteristic spectrum results from the nature of the bipolar explosion. 
For models D5b-D5d (see Fig. 26(b)-(d)), on the other hand, the spectra are quite similar to that of typical type-I 
burst-like gravitational waves [26,29,34]. In these cases, a few sharp peaks appear at the frequency between 1 and 
2 kHz. In Fig. 26(b), the energy power spectrum for smaller-mass model A5b is shown for comparison with the 
spectrum for model D5b. It is found that the shape of the spectra for these two models is quite similar although the 
height of the peak for model A5b is smaller than that for model D5b and the peak slightly shifts to the low frequency 
side due to the fact that the mass is smaller. 

The spectrum for model D5c (see Fig. 26(c)) is similar to that for model D5b except for slight shift of the spectrum 
to low frequency side. The smaller peak frequency results from the fact that the central density at the bounce and 
of the formed protoncutron star for model D5c is smaller than those for model D5b (cf. Figs. 3 and 4), and hence, 
the dynamical time scale is longer. In Fig. 26(c), the energy power spectrum for model D2c is also shown to see 
dependence of the spectra on the angular momentum of the progenitor. The magnitude of the spectrum for model D2c 
is overall smaller than that for model D5c, indicating the magnitude of the angular momentum plays an important 
role for amplifying the gravitational wave amplitude. Also found is that the peak at / w 1 kHz is dominant for model 
D2c. This indicates that for the smaller value of the angular momentum, a single mode is dominantly excited. The 
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FIG. 26. The energy power spectra of A2(f) mode for (a) model D5a, (b) models D5b (solid curve) and A5b (dashed curve), 
(c) models D5c (solid curve) and D2c (dashed curve), and (d) model D5d. Note that gravitational waveforms for models D5a, 
D5b, D5c, and D5d are shown in Figs. 21, 22, 25, and 23, respectively 



reason is likely to be that the degree of deformation of the inner core is smaller for the smaller value of the angular 
momentum, and hence, the fundamental quadrupolc mode of the formed protoneutron star is dominantly excited. 

The spectrum for model D5d (see Fig. 26(d)) has a more complicated shape than those for models D5b and D5c. 
In particular, several small peaks in the low frequency region are found between 100 and 1000 kHz, indicating that 
several oscillation modes arc excited simultaneously. This reflects the manner on the generation of the shocks at the 
surface of the protoneutron star in the equation of state 'd': As mentioned in Sec. IV D, the mass of the formed 
protoneutron star is initially small with the equation of state 'd', and as a result, the accretion of the matter which 
gradually falls from the shock layer subsequently proceeds. The accreted matter intermittently hits the protoneutron 
star and excites its oscillation modes in a complicated manner. As a result, gravitational waves of several characteristic 
oscillation frequencies are generated. 



V. SUMMARY 



We performed fully general relativistic simulations for black hole formation through the collapse of rotating iron 
cores on assumption of the axial symmetry for a wide variety of the mass, the angular momentum, the velocity profile, 
and the equation of state. To systematically study the dependence of the threshold mass for the prompt black hole 
formation on the equations of state, we adopt a parametric equation of state [51,29,34]. We choose the parameters 
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so that the maximum mass of the cold spherical neutron star becomes « 1.6M . The rotating iron cores before 
collapse are simply modeled by the rotating T = 4/3 polytropc in equilibrium. The mass of the iron core is set to be 
« 2.0-3.0M Q . The amount of rotational kinetic energy of the core is changed in a wide range from T Iot /W — to 
« 0.009. 

We have found that the threshold mass for the prompt black hole formation depends sensitively on the angular 
momentum, the rotational velocity profile of the initial core, and the equations of state. The dependence of the 
threshold mass on the these elements can be summarized as follows: (i) The thermal pressure generated by shocks 
increases the threshold mass by 20-40%. The magnitude of this factor depends on the adopted equation of state: (ii) 
With the increase of the spin parameter q, the threshold mass increased by ~ 25g 2 % for the case that the progenitor 
is rigidly rotating: (iii) Effect of differential rotation further increases the threshold mass, since the centrifugal force 
around the central region can be efficiently increased: (iv) The threshold mass depends sensitively on the equations 
of state since the dynamics of the collapse does. For the smaller value of |Ti — 4/3 1 and for the larger value of p n uc, 
prompt black hole formation becomes more liable since the mass of the inner core at the bounce is larger for such 
equations of state. 

About the dynamics of the collapse, we have found the following features: (i) if Ti is close to 4/3, the collapse 
proceeds in an approximately homologous manner, and thus, the most part of the iron core collapses nearly simulta- 
neously. This implies that the mass of the inner core formed at the bounce is larger for the smaller value of |Ti — 4/3 1, 
that helps prompt black hole formation as mentioned above. In the case that the mass of the progenitor is not large 
enough, not a black hole but a protoneutron star is formed. Since the mass of the inner core is larger for the small 
value of |Ti — 4/3 1, the inner core shrinks more significantly resulting in a larger degree of the aspherical deformation 
and in a significant spin-up. As a result, the shape of the shock is highly aspherical and the explosion proceeds in a 
strongly bipolar manner if the progenitor is rapidly rotating as T rot /W « 0.009 (see Sec. IVD 1): (ii) For T x < 1.3, 
the collapse does not proceed in the homologous manner. Instead, only the central region rapidly collapses and 
forms the inner core of a small mass at the bounce. As a result, prompt black hole formation is less liable. For 
the cases of protoneutron star formation, the aspherical deformation of the inner core is not very remarkable since 
its mass at the bounce is not sufficiently large, and hence, the effect of the centrifugal force is not very important. 
Consequently, the shape of the shock is only sightly aspherical even when the progenitor is rapidly rotating (see Sec. 
IVD 1): (iii) For the equation of state 'c' in which the value of p nuc is smaller than that for others, the pressure for 
10 14 g/cm 3 </)<2x 10 15 g/cm 3 is larger than that for the equations of state 'b' and 'd'. As a result, strength of 
the shocks formed at the bounce is enhanced, and therefore, black holes are less liable to be formed. 

Gravitational waveforms are approximately computed in terms of a quadrupolc formula [55]. It is found that the 
amplitude of gravitational waves at the bounce increases monotonically as the mass and angular momentum of the 
iron core increase as far asM< 2.5M . In contrast to the previous results (e.g., [28]), we do not find any tendency 
that the maximum amplitude saturates with the increase of the angular momentum for a given mass. This is due to 
the fact that the general relativistic gravity is strong enough for the inner core to form a compact protoneutron star 
overcoming the centrifugal force in the present choice of the spin parameter (0 < q < 1) for the rigidly rotating case. 

Gravitational waveforms depend sensitively on the equation of state. For models with the equations of state 'b', 'c', 
and 'd', the so-called type-I gravitational waves are emitted even with the mass M ~ 2-2. 5M as in the case of mass 
M w 1.4-1.5-Mq which is studied in the previous papers [34,36]. With the increase of the mass, the amplitude and the 
frequency of gravitational waves become higher. Thus, the difference is only quantitative for the equations of state 'b' 
and V. A point worthy to note is that type-I gravitational waves are emitted for M > 2M even for the equation of 
state 'd' (with Ti = 1.28). For M w 1.5M , type-Ill gravitational waves are emitted for such small value of Ti since 
the mass of the inner core is very small at the bounce. However, in the higher-mass case with M > 2M , the mass of 
the inner core formed at the bounce is large enough to emit type-I gravitational waves. Gravitational waveforms in 
the collapse of rapidly rotating iron core with the equation of state 'a' in which Ti = 1.32 are qualitatively different 
from others. The reason is that in such case, an outstanding bipolar explosion is induced along the rotational axis 
and gravitational waves associated with such extreme explosion and with the resulting oscillation of the inner core in 
the direction of the rotational axis are emitted. Due to this change, the shape of the energy power spectrum is also 
qualitatively different from others. 

In this paper, we assume that the collapse proceeds in the axisymmetric manner. However, the rotating iron core 
may be dynamically unstable against the nonaxisymmetric deformation during the collapse if the spin is increased 
significantly during the collapse, and as a result, the value of T rot /W exceeds a critical value ~ 0.27. To determine the 
criterion for the onset of nonaxisymmetric instabilities, we performed three-dimensional simulation for the iron core 
collapse in a companion paper [38], which shows that in the collapse for rigidly rotating models, the value of T Iot /W 
for the formed protoneutron star is far below the critical value for the onset of the dynamical instabilities, and the 
nonaxisymmetric instability does not set in. Therefore, the assumption of the axial symmetry adopted in this paper 
is justified. 

Finally, we mention a direction of our next study. In this paper, we have clarified a criterion for prompt black hole 
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formation in the iron core collapse adopting the parametric equations of state as a first step. We have shown that 
the dynamics of the collapse and the criterion for the black hole formation depend sensitively on the equations of 
state. To obtain more realistic outputs that will be in nature, in the next step, it is necessary to adopt more realistic 
equations of state. We plan to perform such realistic simulations adopting a realistic equation of state [53] in the fully 
general relativistic framework. 
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APPENDIX A: COMPARISON WITH REALISTIC EQUATION OF STATE 
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FIG. 27. The pressure, P, as a function of the density, p, for the realistic equation of state at zero temperature. The solid 
curve corresponds to the case Y e — 0.5 and the long dashed to Y e = 0.2. The three dotted curves denote the cold equations of 
state with the parameters listed in Table I. 



In this appendix, we compare the parametric equation of state (PEOS) adopted in this paper (see Sec. II B) with 
a realistic equation of state proposed by Shcn ct al. [53] (REOS). In Fig. 27, we show relations between the pressure 
and the density of REOS at zero temperature (adding the degenerate electron pressure) together with those of PEOS. 
Since REOS depends on Y e , we consider the two values Y e — 0.5 and 0.2. The P-p relations of REOS and PEOS 
shows a good agreement, in particular, in the subnuclear density. In the supranuclear density, on the other hand, 
difference between the two equations of state becomes outstanding: The REOS is stiffcr than PEOS (see discussion 
about T 2 below). This fact leads to the larger maximum mass of the spherical cold neutron star for REOS. Indeed, 
the maximum ADM mass is w 2.2M Q [53], much larger than that for PEOS (« 1.6M Q ). 

We first consider the validity of the adiabatic index of PEOS at subnuclear density T\. The adiabatic index Ti 
includes the effect of the electron capture. This implies that it corresponds to an average value (Tm) of the 'effective' 
adiabatic index Tm [85,26] which describes the change of the pressure along a collapse trajectory of a given mass 
element: 
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where 6S = at zero temperature. The crucial quantity for the dynamical behavior of the core is this effective 
adiabatic index Tm [86,26]. 

On the other hand, the P-p relations for REOS in Fig. 27 do not include the effect of the electron capture. Since 
Y e decreases during the collapse due to the electron capture, the P-p relation will shift toward that of the smaller 
Y e . If the electron capture rate is small, the effective adiabatic index will be close to 4/3, while the rate is large, the 
index will deviate from 4/3. 

A plausible value of (Tm) for REOS in subnuclear density would be estimated as follows. First, previous studies 
suggest that the initial value of Y e is w 0.42 and — 0.11 < 8Y e < —0.04 during the collapse [87,13]. Assuming these 
values and 5\ogp « 4, the average value of (5Y e / 5 log p)\m is —0.01 — 0.0275. On the other hand, the value of 
(d\ogP/dY e )\ p , s of REOS is w 1.4-1.9 for 10 8 < p < 10 14 and 0.31 < Y e < 0.42, and its averaged value is « 1.6. 
Thus, the average value for REOS is 1.29 < (r M ) < 1.32 depending the electron capture rate. Therefore, the range 
1.28 < Ti < 1.32 adopted in this paper would be reasonable. 

Then, let us consider the adiabatic index T2 at supranuclear density. The adiabatic index of REOS at supranuclear 
density is larger (r sa 3) than that of PEOS (2.25 < T2 < 2.75). This is due to a relatively larger value of 
incompressibility of the REOS: K s = 281 MeV [53]. Although the real value of K s is uncertain at current status, 
a recent study [88] reported that a plausible range of the value of K s is w 220-270 MeV. Therefore, the adiabatic 
index can be much smaller for a smaller value of K s . For example, the adiabatic index is « 2.2 around the nuclear 
matter density for a realistic equation of state by Lattimer and Swesty with K s = 220 MeV [52]. Thus, the range 
2.25 < T 2 < 2.75 adopted in this paper may be a plausible choice. 

In summary, we conclude that the parameter range adopted in the present paper is not far from predictions and 
suggestion of realistic equations of state. 
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